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The solution of many physical evolution equations can be expressed as an exponential of two 
or more operators acting on initial data. Accurate solutions can be systematically derived by 
decomposing the exponential in a product form. For time-reversible equations, such as the Hamilton 
or the Schrodinger equation, it is immaterial whether or not the decomposition coefficients are 
positive. In fact, most symplectic algorithms for solving classical dynamics contain some negative 
coefficients. For time-irreversible systems, such as the Fokker-Planck equation or the quantum 
statistical propagator, only positive-coefficient decompositions, which respect the time-irreversibility 
of the diffusion kernel, can yield practical algorithms. These positive time steps only, forward 
decompositions, are a highly effective class of factorization algorithms. This work introduce a 
framework for understanding the structure of these algorithms. By a suitable representation of 
the factorization coefficients, we show that specific error terms and order conditions can be solved 
analytically. Using this framework, we can go beyond the Sheng-Suzuki theorem and derive a 
lower bound for the error coefficient evTV ■ By generalizing the framework perturbatively, we can 
further prove that it is not possible to have a sixth order forward algorithm by including only the 
commutator [VTU] = [V, [T, V]]. The pattern of these higher order forward algorithms is that in 
going from the (2n) th to the (2n+2) th order, one must include a new commutator [VT 2n ~ 1 V] in the 
decomposition process. 



I. INTRODUCTION 

Many physical evolution equations, from clas- 
sical mechanics^* 2 * 3 ^, electrodynamics^, statistical 
mechanics^ to quantum mechanics&SiiS, all have the 
form 



dt y 



V)w, 



(1.1) 



where T and V are non-commuting operators. Such an 
equation can be solved iteratively via 



w(t + e) = e t{T+v) w{t), 



(1.2) 
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ovided that one has a suitable approximation for the 

^(T+V) tt„,,„11„ „eT n „A 



short time evolution operator e 



Usually e eT and 



e eV can be solved exactly. By factorizing e e ' T+y ' to 
higher order in the form 



AT+v) 



N 

II' 



,tieT a v z eV 



(1.3) 



one can solve (|1.1[1 accurately with excellent conservation 
properties. Classically, each factorization ll .31) produces a 
symplectic integrator which exactly conserve all Poincare 
invariants. A vast literatur e) 1 ! 2 ! 3 exists on producing sym- 
plectic integrators of the form (|1.3fl . Once a factorization 
scheme is derived, it can be implemented specifically to 
solve any particular evolution equation of the form l|l.l[) • 
However, as one examines these factorization schemes 
more closely, one is immediately struck by the fact that 
beyond second order, all such scheme contain some nega- 
tive coefficients 1 ' 2,3 U and Wj. Since the fundamental dif- 
fusion kernel cannot be simulated or integrated backward 
in time, none of these higher order schemes can be ap- 
plied to time-irreversible systems. This lack of positive- 
coefficient decompositions beyond second order was first 



noted and proved by Shengii. Sheng showed that equa- 
tions for determining the third order coefficients in 1|1.3|) 
are incompatible if the coefficients {ti, Uj} are assumed to 
be positive. This is a valuable demonstration, but it shed 
no light on the cause of this incompatibility nor offered 
clues on how to overcome this deficiency. SuzukiiS later 
proved that the incompatibility can be viewed more geo- 
metrically. His proof tracked the coefficients of the oper- 
ator TTV and TVV in the product expansion of ((Oil . If 
the expansion were correct to third order, then the coef- 
ficients for both operators must be 1/3!. The coefficient 
condition for one corresponds to a hyperplane and the 
other, a hypersphere. Suzuki then went on to show that 
for the same set of positive coefficients, the hyperplane 
cannot intersect the hypersphere and therefore no real 
solution is possible. 

The product form l|1.3[) has the general expansion 

jv , 

JJ e UeT eVi eV = exp I ereT + ev£V + eTy£ 2 y\ 
i=l ^ 

+ e TTV e 3 [T, [T, V]] + e V T V e 3 [V, [T, V]] + ■ 

= e eHA(£ \ (1.4) 

where the last equality defines the approximate Hamil- 
tonian of the product decomposition. The goal of factor- 
ization is to keep ex = ey = 1 and forces all other error 
coefficients such as eyy, cttv, &vtv, etc., to zero. By 
tracing the incompatibility condition to error coefficients 
of specific operators, one can identify which error term 
cannot be made to vanish. The operator TTV can only 
occur in [T, [T, V}} and TVV only in [V, [T, V}}. Thus the 
incompatibility condition is equivalent to the fact that for 
positive coefficients {i;, cttv and eyrv cannot both 
be reduced to zero. To circumvent this, it is suffice to 
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force one error coefficient to zero and keep the other com- 
mutator in the factorization process. Since in quantum 
mechanics [V, [T, V]] corresponds to a local function, just 
like the potential, Suzuki^ suggested that one should fac- 
torize e £ ( T+y ) in terms of T, V and [V, [T, V}}. Following 
up on this suggestion, Suzuki 14 and Chin 1 ^ have derived 
fourth order factorization algorithms with only positive 
coefficients. Chin 1 - also shown that classically, [V, [T, V}] 
give rises to a force gradient exactly as first suggested by 
Ruthi&. Chin and collaborators have since abundantly 
demonstrated the efficiency of these forward time step 
algorithms in solving both time-irreversibloiiiS* 1 *^! and 
time-reversible^iSiiSii^ dynamical problems. Jang et alm^ 
have used these forward factorization schemes in doing 
quantum statistical calculations and Omelyan et alm^^ 
have produced an extensive collection of higher order al- 
gorithms based on this class of fourth order forward al- 
gorithms. 

An important question therefore arises: with the in- 
clusion of the operator [V, [T, V]], can one produce for- 
ward algorithms of sixth or higher order? The answer 
provided by this work is "no" . For a sixth order de- 
composition with positive coefficients, the commutator 
[V, [T, [T, [T, V}]]] cannot be made to vanish and must be 
included. In order to prove this result we have developed 
a formalism to analyze the structure of these forward fac- 
torization schemes. By use of a suitable representation of 
the factorization coefficients, we show that linear order 
conditions and quadratic error terms can both be solved 
analytically. The resulting error term then makes it obvi- 
ous that it cannot vanish if the factorization coefficients 
are purely positive. By use of this formalism we can 
go beyond the Sheng-Suzuki theorem and derive a lower 
bound for the magnitude of the error coefficient eyrv- 
By generalizing the method to sixth order, we further 
prove the main result as stated above. This analytical 
method of solving the order conditions will allows us to 
analyze and classify factorization algorithms in general. 

In the next section we introduce our notations and il- 
lustrate our method of solving the order condition analyt- 
ically by giving a constructive proof of the Sheng-Suzuki 
theorem. In Section III, we discuss the conditions neces- 
sary for a six order forward algorithm. In Section IV we 
introduce a perturbative approach to study the sixth or- 
der case and show that it is not possible to have a forward 
sixth order algorithm by including only the commutator 
[V, [T, V]]. In Section V we discuss the pattern of higher 
order forward algorithms. In Section VI, wc summarize 
our conclusions and suggest directions for future research. 
The Appendix contains details of how to reduce a general 
quadratic error coefficient to a multi-diagonal form. 



II. A CONSTRUCTIVE PROOF OF THE 
SHENG-SUZUKI THEOREM 

In Suzuki's proofs, without explicitly computing Cttv 
and evTVi he showed that both cannot be zero. Here, 



we show that by enforcing exv = and exrv — 0, we 
can compute a lower bound for eyrv analytically and 
show that it cannot vanish for a set of positive {ti}. This 
determination of a lower bound for eyrv goes beyond 
the Sheng-Suzuki theorem in providing a more detailed 
understanding of all fourth order forward algroithms. 

The first step of our approach is to compute the error 
coefficients exvi &ttv e VTV, etc., in terms of the factor- 
ization coefficients {ti,Vi}. This can be done as follow. 
The left hand side of (|1.4|l can be expanded as 

+£ {j2 v ^ v+ ---- ( 2J ) 

Fixing ey = &v = lj the right hand side of l|1.4|) can 
likewise be expanded 

e eH A (e) = i + £ (t + V) + ^e 2 (T + V) 2 + e 2 e TV [T, V] 
+e 3 e V Tv[V, [T, V]] + e 3 e TT v[T, [T, V]] 
+ \e i e TV {(T + V)[T, V] + [T, V](T + V)} 

+ ^e 3 {T + Vf + ■■■. (2.2) 

Matching the first order terms in e gives the primary 
constraints 

N N 

^t i = l and ^Ui = l. (2.3) 

i=l i=l 

To determine the other error coefficients, we focus on 
a particular operator in l|2.2|l whose coefficient contains 
ctv, &ttv or eyrv and match that operator's coeffi- 
cients in the expansion of l|2.1|) . For example, in the 
e 2 terms of (|2.2|l . the coefficient of the operator TV is 
(i + exv)- Equating this to the coefficients of TV from 
(|2T|) gives 

1 N 

- +e T y =y]sjVi, (2.4) 
1 i=i 

where we have introduced the variable 

i 

Si = J2tj- ( 2 -5) 

Alternatively, the same coefficient can also be expressed 
as 

1 N 

77 + e VT = Uuj. (2.6) 

i=i 

where 

N 

U i = ^2 v r ( 2J ) 
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It turns out that Sj and U{ are our fundamental variables, 
the coefficients tj and x>i are backward and forward finite 
differences of s» and Ui, 



U = Si - Sj_i = Vsj 
= Ui -u i+ i = -Viii 



(2. 



The results l|2.4l) and (|2.6(l are equivalent by virtue of the 
"partial summation" identity 



N N 

2J Vs l u l = - 2J SiAui. 

i=l i=l 



(2.9) 



(Note that sq = and un+i = 0.) In the following, 
we will use the backward finite difference operator exten- 
sively, 



Vs? 



(2.10) 



with property 



N 

E 

i=l 



Matching the coefficients of operators TTV and TVK 
gives 



1 1 



N 



N 



+ -<t\ +e TTV = = -^Vs^i, (2.11) 

i=l i=l 
1 * 

= -E Vs ^- ( 2 - 12 ) 



1 1 

— I + 2 e TV - ervT 2 



i=i 



The error coefficient eyrv can be tracked directly by the 
operator VTV. The coefficient for the operator VTV is 
quadratic in Vi but not diagonal. This is more difficult 
to deal with than TVV's coefficient. Nevertheless, we 
show in the Appendix that, VTV's coefficient can be 
diagonalize by a systematic procedure to yield the same 
constraint equation as H2.12fl . 

In order to have a fourth order algorithm, aside from 
the primary constraints (|2.3[1 . one must require exv = 0, 
&ttv — 0, and eyrv — 0. F° r a symmetric product form 
such that t\ = and Vi = Ujv-i+i, ij+i = ijv-i+i) ° r 
vn = and Vi — v^-i, ti — tw-i+i, one has 



c -eH A (~e) e eH A (e) _ ^ 



(2.13) 



This implies that Ha(s) must be a even function of e, 
and exv = is automatic. The vanishing of all odd 
order errors in Ha{£) implies that we must have 



(2n- 1) 



i 



{2n)\ 



(2.14) 



ensuring that T 2n ~ x V has the correct expansion coeffi- 
cient. It is cumbersome to deal with symmetric coeffi- 
cients directly, it is much easier to use the general form 



(|1.3|) and to invoke i|2.14|) when symmetric factorization 
is assumed. 

The next step in our strategy is compute a lower bound 
for the magnitude of eyTV ', after satisfying constraints 
ctv — and errv = 0. We view latter two constraints 



i=l 

N 1 
E VS ^ = Q ' 



(2.15) 
(2.16) 



as constraints on {ui} for given a set of {ti} coefficients. 
For positive {ti}, the RHS of l|2.12|l is a positive-definite 
quadratic form in Wj. Its lower bound can be deter- 
mined by the method of constrained minimization using 
Lagrange multipliers. Minimizing 

1 N ( N l\ 

f = ^y, Vs ^ ~ A i E Vs ^~2 

i=\ \i=l / 



gives 



Vs ■ Vs 2 

Ui = Ai— - + A 2 — i = Ai + A 2 (s, + Si_i). (2.18) 

VS; VS; 



Imposing (|2.15|) and l|2.16|l determines Ai and A 2 , 



Ai + A 2 — -, 



Ai + A 2 + g\ 2 = -, 



where g defined by 



^VsfVs 2 



Vs ? ; 



= 1+3, 



is given by 



A' 



= ^ SiSi-l(si - Si_l). 



(2.19) 
(2.20) 

(2.21) 
(2.22) 



By substituting in SjS;_i = (s? + sf_ 1 - (s, - Si_i) 2 )/2, 
one discovers that 

.9 = ~2 g+ 2^ ~ ^ 

and therefore 

1 - 
<? = -(!-£<?), where fc = £*i- (2.23) 



The factor 1/3 is the continuum limit (N — > cx)) of g 
when the sum is replaced by the integral f Q s 2 ds. The 
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evaluation of general sums of the form i|2.22|l will be fur- 
ther discuss below. This exact form for g obviated the 
need to determine g's upper bound as it is done origi- 
nally in the work of SuzukiiS, and in the more recent 
work on symplectic correctors*.) With Ai and A2 known, 
the minimium of F is given by 



F = 



(Ai + A 2 ) 2 + -g\\ 



and therefore, 



1 

2 

1 1 

4 + T2g 



evTV ^ 



69 



24(1-^)' 



(2.24) 



1 



Sg 



24(l-6g) 



(2.25) 



This implies that, first, eyrv must be negative, 
ondly, its magnitude is 



Sec- 



\evTv\ > 



1 



Sg 



24(l-5g) 



(2.26) 



The Sheng-Suzuki theorem now follows as a simple corol- 
lary. If all the Vs are positive, then eyrv cannot vanish 
because its lower bound H2.26J) . which depends on 5g as 
given by l|2.23[l . cannot vanish. The only way to achieve 
a fourth order forward algorithm is to keep the commu- 
tator [V, [T, V]] with coefficient evrv, but move it to the 
left hand side of i|1.4|) . This means that for all such fourth 
order algorithms, the sum of factorization coefficients of 
all the [V, [T, V]] terms must be positive. All such fourth 
order algorithms are characterized by their respective val- 
ues of evrv 5 and how well they saturate the lower bound 
(|2.26|) . Note that in deriving this lower bound, we did 
not need to incorporate the primary constraints u\ = 1. 

A very different "elementary" proof of the Sheng- 
Suzuki result has been offered by Blanes and Casa^i 
Our work is more precise in demonstrating that, not only 
evrv cannot vanish, it has a lower bound (|2.26|) deter- 
mined only by {ti}. 

Note also that Vi = ui — Ui + i and (|2.18fl implies that 



Vi = A 2 (si_i 



1 fa + t i+1 ) 

2 (1-Sg) ■ 



(2.27) 



Thus, if one insists that eyTV be zero, then Sg can be zero 
only if at least one ti is negative such that (ti + or 
(tj +i;_i) remains negative. Ea. (|2.27|) then implies that 
its adjacent values of Vi or vi-\ must also be negative. 
Thus a fourth order factorization without keeping any 
additional operator such as [V, [T, V]] must have at least 
one pair of negative U,Vi coefficients. This result was 
first proved by Goldman and Kaper— . This simpler proof 
follows the idea of Blanes and CasaSi 



III. THE SIXTH ORDER CASE 

By incorporating the potential-like operator [V, [T, V]], 
many familiesi2*22i2£ of fourth order forward algorithms 



have been found. They are not only indispensable for 
solving time-irreversible problemsi£ii2ii2i2S; they are also 
superior to existing fourth order algorithms in solving 
time-reversible classical^i^SSiS and quantum2i±& dynam- 
ical problems. It is therefore of great interest to deter- 
mine whether there exist practical forward algorithms of 
even higher order. We show in this section that sixth or- 
der forward algorithms requires the inclusion of the com- 
mutator [V, [T, [T, [T, V}]]]. The inclusion of [V, [T, V]] 
which make possible fourth order forward algorithms, is 
insufficient to guarantee a sixth order forward algorithm. 
In general, if F2„(e) is a 2nth order forward decompo- 
sition of e £ ( T+v \ then F 2n +2(£) would require the in- 
clusion of a new operator not previously included in the 
construction of F2 n (e). We have proved the case of n = 1 
in the last section. The new operator being 



V 1 = [V,[T,V]}. 



(3.1) 



Consider now the case n — 2. In the following dis- 
cussion, we will use the condensed bracket notation: 
[V 2 T 3 V] = [V, [V, [T, [T, [T, V}}}}}, etc.. We have shown 
in the last section that, for positive ti, with Ui satisfying 
constraints (|2.15l) and l|2.16|l . we can factorize e £ ^ T+v ^> 
up to the form 



N 

n 



e UeT e v iS V 



exp 



sir + V + e VTV e 2 [VTV] 

4 -, 

+e 4 ^ ei Q 4 + 0(£ 6 )) , (3.2) 



where eyrv cannot be made to vanish, and Qi are 
four independent operators described below. There is 
one error operator [TV] in first order, two error opera- 
tors [TTV] and [V7V] in second order, four operators 
[TTTV], [VTTV], [TVTV] and [VVTV] in third order, 
and eight operators 

[TTTTV], [VTTTV], [TVTTV], [VVTTV], 
[TTVTV] , [VTVTV], [TVVTV] , [VVVTV], 

in fourth order. These error operators are results of con- 
catenating T and V with lower order operators on the 
left. In each order, not all the operators are independent. 
For example, setting C = [AB] in the Jacobi identity 

[ABC] + [BCA] + [CAB] = 0, 

gives [ABC] = [BAG] and therefore 

[ABAB] = [BAAB]. 

For the case where [FTV^] commutes with V we also have 
[]/ n ]/TT/] = 0. Hence there are only two independent 
operators [TTTV], [TVTV] in third order and four op- 
erators [TTTTV], [VTTTV], [TTVTV] , [VTVTV] in 
fourth order. The last two are just [TTVi] and [VTVi], 
which resemble second order errors for a new potential 
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V\. To have a sixth order algorithm, one must eliminate 
these four error terms. Since [TTVi] and [VTVi] are 
linear in Vi, they can always be eliminated by including 
sufficient number of V\ operators in the factorization pro- 
cess. The remaining error terms [T 4 T^] and [VT 3 V] are 
unaffected by V\ and can only be eliminated by the choice 
of coefficients {ti,Vi}. Thus we can apply our previous 
strategy of dealing only with coefficients {U, v{\ but now 
computing the error coefficient evr 3 v explicitly. 

A careful reexamination of our proof for the Sheng- 
Suzuki theorem shows that we have proved more than 
that's required. The minimization procedure produces a 
lower bound for evrv, whereas the Sheng-Suzuki theo- 
rem only requires that evrv not be zero. The expansion 
()2.18(l merely served as a vehicle for demonstrating that, 
for any {ui} satisfying l|2.15[l and i|2.16[) . evrv cannot 
vanish for positive {ti}. We do not really need to min- 
imize anything, or to determine an actual lower bound. 
This suggests a simple strategy for proving the sixth or- 
der case. It is sufficient to show that eyT 3 v cannot vanish 
for any set of {u{\ satisfying higher order constraints. 



= 1 + 9n 

where we have used the identify 



Y7s m Y7« n V7s m_1 Vq n_1 
VS i VS i _ V7„™+"-l ' N ' 



Vs,; 



SiSi-i- 



(4.5) 



(4.6) 



to define the reduced symmetric matrix g mn - Since 
Gin = G n \ = 1 (and hence g\ n — g n \ = 0), we can 
subtract the first constraint equation 



Ai + A 2 + A3 + A4 



1 



(4.7) 



from the other three and reduce the system down to three 
equations for m — 2 to 4: 



V - 1 1 

m + 1 2 

n— 2 



(4.8) 



By writing, s 2 = s,_i + ±Vsj and Sj_i = s 4 _i - \Vsi 
where s i _i = |(s, + s,-_i), we can systematically expand 



IV. PROVING THE SIXTH ORDER CASE 

As discussed in the last section, for a sixth order al- 
gorithm, a symmetric factorization must satisfy, in ad- 
ditional to (|2.15|l and (|2.16f) . the constraint I|2.f4|l for 

n = 2, 



1 



N 

»=i 



(4.1) 



Also, since the operator T V uniquely tracks the commu- 
tator [T 4 !/], the error coefficient e T i V will vanish if the 
expansion coefficient of T 4 V is 1/5!. This means that 
factorization coefficients {U,Vi} must also obey 



1 



N 

<=1 b 



(4.2) 



These four constraints (|2~T51) . (|2~To|) . JOJ, and (|4~2l . can 

be satisfied by the expansion, 

x Vs 2 , Vsf , Vs 4 
Ui = A 1 + A 2 — ± + \ 3 —± + \ 4 —±. (4.3) 

VSi VSi VSi 

We must now demonstrate that in this case, eyT a v can- 
not vanish if {ti} are all positive. 

When Ui is expanded via (|4.3I) . the four constraints 
(|2~T5|) . (|2~To|) . JO}, and (|4^|) produce the following set 
of four linear equations for m — 1 to 4, 



7 , G mn X n . 

m + 1 

n— 1 



(4.4) 



The matrix G mn is given by 

N 



N 



Gmn — ^ — — 1 + SiSi-i 



i=l 



i=l 



l!(n-l)! 



3!(n-3)!2 2 »'-* 



When each summant Vs™Vs"/Vsi is expanded and com- 
pared with the similarly expanded integral 



s m+n-2 ds = 



' Si-l 

we deduce that 

G r nn. 



m + n — 1 



-Vs 



m+n— 1 



m + n — 1 



1 f w 

--mn(m - l)(n - 1)^ $>£$"- 4 (V*0 8 
I »=i 

AT 

+ A 5 ^s™1"- 6 (Vs l ) 5 + --- 



,(4.9) 



with 



A = — [(m+n-4) 2 + (m-2)(2m-7) + (n-2)(2n-7)]. 

The constant part of the matrix is the continuum limit 
(JV — > 00) of the sum, which is the integral 



mns 



m+n— 2 



ds 



mi) 



m + n — 1 



We will denote this constant part of the matrix as G mn . 
The corresponding continuum part of g mn is g mn = 
G mn ~ 1- The remaining finite parts of G mn in (|4.9|l . 
which depends explicitly on Sj, will be denoted as SG mn - 
Since g mn differs from G mn only by a constant, its finite 
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part Sg mn is the same as that of G mn , i.< 
By repeated applications of the identity 
duce g mn to a sum of terms of the form 



one can re- 



K(l 7 n) = J2( s ^y Vs i- 



(4.10) 



Since the explicit form of g mn is known via 14. 9|) , these 
functions are not particularly useful as calculational 
tools. However, they are very useful in quickly identi- 
fying the matrix element of g mn when doing analytical 
calculations. For later reference, we list below some 5 mn 's 
in terms of k(1, n): 



922 = 


k(1,1) 




923 = 


k(1,2) 




524 = 


k(1,3) 




.932 = 


k(1,3)H 


-k(2,1) 


533 = 


k(1,4)H 


-k(2,2) 


534 = 


k(1,5)H 


-k(2,3) + «(1,3) 



(4.11) 



Note that 1722 is the g function of the last section. From 
the general formula (|4.9|) , one finds indeed that g 22 = 1/3 
and 



692 



1 N 



1 



(4.12) 



1=1 



If we only keep the continuum matrix g® nn in 14. 8|) 



the solution is trivial: A2 




1 

6 
_L 

~4 
J_ 
10 



\ 2 = — J, -^3 = 0, A4 = 0. This 
suggests that we should also expand each A, into its 
continuum and finite part: A2 = — \ + SX 2 , A3 = <5A 3 , 
A4 = 6X4. For our purpose, it is enough to keep the lead- 
ing finite size correction term, i.e., we can neglect the 
terms of the form Sg mn SXk- In this case, we just have 



0, X 4 




sx 2 \ 


1 \5g 2 2 


8X 3 = 


5*523 


5Xj 


\ 2-^524 



(4.13) 



We do not need to solve each SXk explicitly; we only need 
to know that they are proportional to Sg 2n . Since Ai + 
A2 + A3 + A4 = |, this also implies that Ai = 1 + <5Ai with 



SXx + 5X 2 + SX 3 + 5X 4 = 0. 



(4.14) 



The above discussion suggests that one should also sep- 
arate Ui into its continuum and finite part, 



(1 



l_Vs| 

2 Vs, 



) + Sui 



(4.15) 



The constraints on u, now translate into constraints on 
(5m, : 



N N 
i=l i=l 



2 Vs,' 

1 



T - (1 - ^G 2n ) = -Sg 2n (4.l6) 



Recall that since g\ n — g n \ — 0, we also have Sg n i — 
Sgin = 0. The above constraints for Sui is exact. We 
have not yet invoked any particular representation for 
5ui. 

To illustrate how this formalism will be used, let's re- 
compute the quadratic form of the last section: 



JV 



X>i«?=E v * (1--^)+^ 



1=1 »=i 

N 



1V^| 
2 Vs,- 



JY 



N 



1 V« 2 Y7« 2 

+ 2 ]T V* 6^ - ]T Vsl^i + 0(6i 



4 f-f Vs, 

1=1 4 — 1 2 — 1 



1^ L 1 1. 11 
i G 22 - -Sg 22 = - - -^22 = 3 + ^5- 



This then implies that 



1 



N 

24^*i 

i=l 



(4.17) 
(4.18) 



(4.19) 



The first key observation is Eq. 1)4. 17|) : to leading order in 
3g2n, this quadratic form only depends on the first two 
constraints on Suj,. Its leading finite part is unchanged by 
additional, higher order constraints on Sui. That is, Sui 
can be very general. By inspection, eyrv above cannot 
vanish for positive {£,}. Thus this leading order calcu- 
lation, while not sufficient to determine the exact lower 
bound for evTV, it is sufficient to show that evTV cannot 
vanish, and thus proving the Sheng-Suzuki theorem. 
Secondly, if Sut were to be represented as 

Vs 2 Vs 3 Vs 4 

Su i =SX 2 (^-l)+SX 3 (^-l)+SX 4 (^-l), (4.20) 

VSi VSi VSi 

then in order for the constraints l|4.16[) to determine SXk 
to the same leading order in Sg 2n as in (|4.13l) it is enough 
to compute only the constant (continuum) part of any 
sums multiplying SXk- This implies that we may replace 
any such sum by its integral, or by any other sum having 
the same integral. Thus for any sum multiplying Sui, we 
may replace it by another sum having the same integral. 
This crucial simplification makes it unnecessary to solve 
for each Xk explicitly. 

To compute the error coefficient eyT 3 v, one must 
use an operator that tracks the commutator [VT 3 F] 
uniquely. The analogous operator T 3 V 2 , whose expan- 
sion coefficient is easy to compute, is no longer suitable. 
Let Ctzv 2 denote its expansion coefficient in terms of 
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{U,Vi} from the left-hand-side of l|3.2|) . By matching 
the same operator's expansion coefficient from the right- 
hand-side, one finds 2 -* 

C T 3 V 2 = — — y e VTV ~ &T 2 VTV - (4-21) 

It is difficult to disentangle evT 3 v from the contaminat- 
ing effects of evTV and e^vTV- The three operators that 
track [VT 3 V] uniquely are VT 3 V, VT 2 VT, and TVT 2 V . 
We choose the symmetric choice VT 3 V, whose coefficient 
is related to eyr^v by 



C VT 3 V = — + 2e VT 3 V . 



(4.22) 



From the left hand side of l|3.2[l . one deduces 

N-l N 
i— 1 j— z+1 

This quadratic form in {^i} is difficult to work with be- 
cause it is not diagonal in Ui or some other variables. In 
the Appendix, we show that it can be simplified to the 
following bi-diagonal form, 

1 / N N i\ 

' \ z=l i=l / 



where z% is defined by 



= £«. 



(4.25) 



Given the expansion l|4.3() for u,, we can deduce the 
corresponding expansion for Zj. From (|4.25|) . we can 
rewrite as 



Zi = U. L Si + U j^ S 'j 
3 =l+l 



(4.29) 



For Ui — A„Vs"/Vsi, we have 
"Vs? 



= Ar, 

= Ar. 



Vs.. 



■Si + (1 - *?) 



1 + SjSi-i 



V.s 

vsr 1 - 



> + (!-*?) 



(4.30) 



Hence corresponding to (|4.3I) , Zj has the expansion 

Vs 2 

= Ai + A 2 (l + s i Si_i) + A 3 (l + SjS i _i— -) 

Vsi 

+\ 4 (l + Si s^A). (4.31) 

One can check that this form for Zi satisfies the four 
constraints 12.15|l . I|2.16|l . (|4.1() . and (|4.2|l when they are 
expressed in terms of Zf 



z\ — Ai + A 2 + A3 + A 4 



1 

2' 



and for m = 1 to 3, 



E v -f 



1 



m + 2 



(4.32) 



The required coefficient ei/T 3 v can now be computed 
from 



e VT 3 V 



/ N N 
\ i=l i=l 



fu 2 -- . (4.26) 



The quadratic form involving u\ is 

^ V*?«? = J2 Vaf (1 - i |^) 2 + 2 ^ Vs? *u 



(=i 



JV 



+ 0(*u?) (4.27) 



3 1 3 

j - G 32 + ^(G 33 + G 24 ) + Sg 2 3 - -^Sg 2 4 



1 l r l r 
10 + 4^ 33 " 2 Sm - 



(4.28) 



In 14.27|l . we have replaced the sum involving Vsf Vs 2 /Vs, 
by its integral equivalent (3/2)Vsf. Also, we have used 
the identity 



Vsf /Vs 2 Vsf\ _ Vsf Vsf + Vs|Vs 
Vsi I Vsi / Vsj Vsi 



Vs.? 



The identity (|4.6|) is needed to show that (|4.32|) is equiv- 
alent to the last three constraint equations for ttj. As in 



the case of ut, we can write z^ in the form 



Zi = -(1 - SiSi^) + SZi 



(4.33) 



and transfer the last three constraints on Zi to 5zi, 

N 



E Vs " 1<5z * = ^52«- 



(4.34) 



The quadratic form for Zj is then 

AT ^ JV 

4 



VsjZ 2 = - X Vsj(l - s 4 s 4 _i) 2 + ^2 Vsi^i 

z— 1 i— 1 i— 1 

iV 

-^2 V8i(siSi-i)5zi + O(Szf) 

i=l 

~ - 1) + i«(2, l) + i<5 322 - i ^ Vsf ^ 

~ 2 — 1 

11 1, , 1. 1. 

I - 2-922 + I (.933 - 524) + 2-6322 - g0ff24 

^ + J-5.933 - ^.924 (4.35) 
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We have again replaced the sum involving Vsi(siSi_i) by 
its integral equivalent (l/3)Vs 3 and used l|4.11|l to express 
the required sum in terms of g mn 's. Thus the bi-diagonal 
form is 

N N 

3^ Vs^ 2 - V*Ni = — + 7 (25g 33 - 36g 2i ) 



i=i »=i 
From (|4.9|l we find, 

N 



10 4 



i=l i=l 
JV jV 

5.924 = -2^^_i( Vs ') 3 - i^E( Vs *) 5 > ( 4 - 36 ) 



i=l 

and therefore finally, 



i=l 



JV iV 

^=24oB^ = ^E^ (4-37) 

i=l z=l 

This is remarkably similar to l|4.19|l . Thus if {ti} are all 
positive, then evT 3 v cannot vanish. No sixth order pos- 
itive factorization scheme is possible without including 
the commutator V 3 = [VT 3 V]. 



V. BEYOND SIXTH ORDER 

In Sections II, we have shown that in order to have 
a fourth order forward algorithm, one must include the 
commutator Vi — [VTV] in the factorization process. In 
the last section, we have proved that in order to have a 
sixth order forward algorithm one must include in addi- 
tion to Vi, the commutator V 3 = [VT 3 V]. By repeating 
the same argument, it is not difficult to discern the pat- 
tern of higher order forward algorithms. In going from 
the (2n) th to the (2n+2) th order, one must add a new 
commutator 

y 2 „_! = [VT 2n - l V] 

to the factorization process. A proof of this general result 
is a straightforward generalization of our approach in the 
last section, but technically much more involved. For 
example, to prove the eighth order case, we must track 
e V"r 5 v uniquely via the operator VT 5 V's coefficient given 
by 55/5!, where S5 as shown in the Appendix, is tri- 
diagonal in m, Zi and 



Vi 



N 
3=i 



One then has to work out the expansion for y, as in the 
case of Zi. Moreover, since eyT 5 v is anticipated to be 
oc J2iLi(^ s i) 7 1 one can n0 longer ignore contribution of 
order (Sui) 2 oc {J2iLiC^ s i) 3 ) 2 ■ Thus the current formal- 
ism, while powerful in determining eyrv variationally 
and evT 3 v perturbatively, is too demanding for the gen- 
eral case. To prove such a general result, one must find 
a less explicit approach. 



VI. CONCLUSIONS 



In this work, we have introduced a framework for an- 
alyzing and understanding the structure of factorized al- 
gorithms. There are three key ideas: 1) The order con- 
straints and error coefficients can be tracked by operators 
and expressed directly in terms of factorization coeffi- 
cients. 2) By introducing a suitable representation for 
the factorization coefficients, the order constraints and 
error terms can be solved analytically. 3) For many pur- 
poses, it is sufficient to determine the error coefficients 
perturbatively. This last point is specially important. All 
previous works on factorization algorithms are based on 
exact decompositions. Since this is difficult to do analyt- 
ically, one can make little progress except numerically . 
This work shows that a leading order calculation is suffi- 
cient to establish most of the important results we know 
about these algorithms. In particular, we have provided 
a constructive proof of the Shang-Suzuki theorem. Most 
importantly, we have shown that in order to have a sixth 
order forward time step algorithm, one must include the 
commutator [yT 3 y] in the factorization process. 

This work suggests that there is regularity to the exis- 
tence of forward algorithms. In order to have only pos- 
itive time steps, one must continue to enlarge one's col- 
lection of constituent operators for factorizing e etyT+v \ 
For a (2n) th order forward algorithm one must include 
all commutators of the form [l/T 2fc_1 V / ] from k = 1 to 
k = n — 1, in addition to T and V. The proof of this 
general result is currently beyond scope of our perturba- 
tive approach. Moreover, the massive cancellations that 
produced the sixth order result l|4.37[) strongly suggest 
that a better formulation, with these cancellations built- 
in, must be possible. This work suggests that a more 
powerful way of understanding the structure of these al- 
gorithms is still waiting to be found. 

The need to include [FT 3 !/] make it difficult to con- 
struct, but does not necessarily preclude the possibility of 
a sixth order forward algorithm. One simply has to work 
harder to devise practical ways of obtaining [VT 3 F] with- 
out computing it directly. Work is currently in progress 
toward this goal. 
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APPENDIX A: COEFFICIENT OF VTV, VT 3 V 
AND VT 5 V 

There is a systematic way of diagonalizing the sum 

N-l N 

Sm = E E Vi ( s i ~ Si ) mv 3 

i—l j=i-\-l 

needed in computing the error coefficients evT m v- The 
above is a sum over the upper triangle of a N x N square 
matrix and can be denoted more simply as X^>«- 
The general form we need to diagonalize is 



complicated operator VTV determines the same evTVi 
as it must. 
For m = 3, we have 

S 3 = E^( s j - s%)vj - 3^2viSi(sj - SijsjVj 



3>i 



j>i 



Assuming now that all linear constraints on m are sat- 
isfied up to the relevant order, we have for the first and 
second term on the right respectively, fi = Vi, gt = sf, 
hi = m, F = 1, P = \ and fi = SjWi, gi = Si, hi = Zi, 
F = g, aP = -|. Hence we have 



JV N 

S(f, g) = J2 Mi ~9i)fj = E fwfi-Y,M ^ A1 ) S " = 4 " V ^ ~ 3 (e " ? Vs 



z 2 \ 



l>3 



j>i 



where we have interchanged the summation indices in the 
first term on the right-hand-side. The key point here is 
that if we introduce a new variable 



N 



hi = ^2fj, 



3=i 



such that fi = hi — /ij+i, then the second term on the 
right hand side of IjAlfl is only a single sum. The first 
term can be eliminated by completing the "square ma- 
trix" . Let Yli fidi — P an d J2j fj — F be known sums, 
then 

pp = J2 E fi = E fa + E fi9ifi + E foiifi- 

i j i i>j j>i 

(A2) 

Subtracting (|A1I) from l|A2[) gives 

PF - S(f, q) = E fhi + 2 E faifi 



N 



N 



Efi(^-i - h l+ i) 2 + 2^2gi(hi - h l+1 )h l+1 

i=l i=l 
N N 

E^-^+i)=E V ^A 2 I (A3) 



and hence, 



A ! 



S(f,g) = PF-J2^9ih 2 



(A4) 



For the case of m — 1, we have fi = m, gi = Sj, hi — 
m, F = 1 from JOl, and P = (i + e T y) from itOjl . 
Therefore, we have 



1 " 
5i = (- + e T v) - E Vs »"i 

i=l 



ul 



Since the coefficient of VTV is just = ^ + e^y, 
the above is identical to (|2.12|) . The use of the more 



where 



AT 

Zi = V <'•;> 



The coefficient of VT 3 V is 5 3 /3!. Since [VT 3 y] contains 
the operator VT 3 V twice, we have 

^^3 = 77 + 2ey T 3y, 

6 5! 



and therefore 

12e VT3v = 5 3 -7^ = 3E V Sl ^-E Vs ' u ?-^ ( A5 ) 



iV 



i=l i=l 

For the case to = 5, we have 



j>i j>i 

+10E«i4(sj - (A6) 

For the first term we have fi = Vi,gi = sf,hi=Ui,F = 
1, and P = j|. For the second term we have /j = Sj«j, 
<7i = s| , hi — Zi, F — and P = g. For the third term, 
we have fi = sf«i, <fc = Sj, ft.^ = y 2 , F = |, and P = \- 
We therefore have 

^ = ^-E v ^- 5 (^-E Vs ^ 2 

2 — 1 2 — 1 

1 ^ 

+ 10 (t2"£ VS4; ' 

i=l 

Af AT AT 

= i-E Vs ^ 2 + 5 E Vs ^- 10 E Vs ' 



i=l 



where 



;.=i 



N 

Vi = J2 v i s l 

3=i 



i=l 
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